sink("Log_Rep_RDD.txt")

# Lampedusa 1

data61$lamp <- as.numeric(data61$date) - 15981
hist(data61$lamp, breaks = 1000)
rdlamp1 <- RDestimate(factor ~ lamp, data = data61, cutpoint = 0, bw = 30)
summary(rdlamp1) 

# with region fixed effects
summary(RDestimate(factor ~ lamp | region, data = data61, cutpoint = 0, bw = 30))

# with covariates
covlamp1 <- RDestimate(factor ~ lamp | education + gndr + agea + foreignparent, data = data61, cutpoint = 0, bw = 30)
summary(covlamp1)


data62$lamp <- as.numeric(data62$date) - 15981
hist(data62$lamp, breaks = 1000)
rdlamp2 <- RDestimate(factor ~ lamp, data = data62, cutpoint = 0, bw = 30)
summary(rdlamp2)

# with region fixed effects
summary(RDestimate(factor ~ lamp | region, data = data62, cutpoint = 0, bw = 30))

# with covariates
covlamp2 <- RDestimate(factor ~ lamp | education + gndr + agea + foreignparent, data = data62, cutpoint = 0, bw = 30)
summary(covlamp2)

low61 <- rdlamp1$ci[1,1]
up61 <- rdlamp1$ci[1,2]
coef61 <- rdlamp1$est[1]
se61 <- rdlamp1$se[1]

low62 <- rdlamp2$ci[1,1]
up62 <- rdlamp2$ci[1,2]
coef62 <- rdlamp2$est[1]
se62 <- rdlamp2$se[1]

# 31.08.2014
data71$one <- as.numeric(data71$date) - 16313
hist(data71$one, breaks = 1000)
rd1 <- RDestimate(factor ~ one | cntry, data = data71, cutpoint = 0, bw = 30)
summary(rd1)

low1 <- rd1$ci[1,1]
up1 <- rd1$ci[1,2]
coef1 <- rd1$est[1]
se1 <- rd1$se[1]

# with region fixed effects
summary(RDestimate(factor ~ one | region, data = data71, cutpoint = 0, bw = 30))

# with covariates
cov1 <- RDestimate(factor ~ one | cntry + education + gndr + agea + foreignparent, data = data71, cutpoint = 0, bw = 30)
summary(cov1)

data72$one <- as.numeric(data72$date) - 16313
hist(data72$one, breaks = 1000)
rd1b <- RDestimate(factor ~ one | cntry, data = data72, cutpoint = 0, bw = 30)
summary(rd1b)

low1b <- rd1b$ci[1,1]
up1b <- rd1b$ci[1,2]
coef1b <- rd1b$est[1]
se1b <- rd1b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ one | region, data = data72, cutpoint = 0, bw = 30))

# with covariates
cov1b <- RDestimate(factor ~ one | cntry + education + gndr + agea + foreignparent, data = data72, cutpoint = 0, bw = 30)
summary(cov1b)


# 14.09.2014
data71$two <- as.numeric(data71$date) - 16327
hist(data71$two, breaks = 1000)
rd2 <- RDestimate(factor ~ two | cntry, data = data71, cutpoint = 0, bw = 30)
summary(rd2)

low2 <- rd2$ci[1,1]
up2 <- rd2$ci[1,2]
coef2 <- rd2$est[1]
se2 <- rd2$se[1]

# with region fixed effects
#summary(RDestimate(factor ~ two | region, data = data71, cutpoint = 0, bw = 30)) # NaNs produced 

# with covariates
cov2 <- RDestimate(factor ~ two | cntry + education + gndr + agea + foreignparent, data = data71, cutpoint = 0, bw = 30)
summary(cov2)

data72$two <- as.numeric(data72$date) - 16327
hist(data72$two, breaks = 1000)
rd2b <- RDestimate(factor ~ two | cntry, data = data72, cutpoint = 0, bw = 30)
summary(rd2b)

low2b <- rd2b$ci[1,1]
up2b <- rd2b$ci[1,2]
coef2b <- rd2b$est[1]
se2b <- rd2b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ two | region, data = data72, cutpoint = 0, bw = 30))

# with covariates
cov2b <- RDestimate(factor ~ two | cntry + education + gndr + agea + foreignparent, data = data72, cutpoint = 0, bw = 30)
summary(cov2b)

# 08.02.2015
data71$three <- as.numeric(data71$date) - 16474
hist(data71$three, breaks = 1000)
rd3 <- RDestimate(factor ~ three | cntry, data = data71, cutpoint = 0, bw = 30)
summary(rd3)

low3 <- rd3$ci[1,1]
up3 <- rd3$ci[1,2]
coef3 <- rd3$est[1]
se3 <- rd3$se[1]

# with region fixed effects
summary(RDestimate(factor ~ three | region, data = data71, cutpoint = 0, bw = 30))

# with covariates
cov3 <- RDestimate(factor ~ three | cntry + education + gndr + agea + foreignparent, data = data71, cutpoint = 0, bw = 30)
summary(cov3)

data72$three <- as.numeric(data72$date) - 16474
hist(data72$three, breaks = 1000)
rd3b <- RDestimate(factor ~ three | cntry, data = data72, cutpoint = 0, bw = 30)
summary(rd3b)

low3b <- rd3b$ci[1,1]
up3b <- rd3b$ci[1,2]
coef3b <- rd3b$est[1]
se3b <- rd3b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ three | region, data = data72, cutpoint = 0, bw = 30))

# with covariates
cov3b <- RDestimate(factor ~ three | cntry + education + gndr + agea + foreignparent, data = data72, cutpoint = 0, bw = 30)
summary(cov3b)

# 18.04.2015
data71$four <- as.numeric(data71$date) - 16543
hist(data71$four, breaks = 1000)
rd4 <- RDestimate(factor ~ four | cntry, data = data71, cutpoint = 0, bw = 30)
summary(rd4)

low4 <- rd4$ci[1,1]
up4 <- rd4$ci[1,2]
coef4 <- rd4$est[1]
se4 <- rd4$se[1]

# with region fixed effects
summary(RDestimate(factor ~ four | region, data = data71, cutpoint = 0, bw = 30))

# with covariates
cov4 <- RDestimate(factor ~ four | cntry + education + gndr + agea + foreignparent, data = data71, cutpoint = 0, bw = 30)
summary(cov4)

data72$four <- as.numeric(data72$date) - 16543
hist(data72$four, breaks = 1000)
rd4b <- RDestimate(factor ~ four | cntry, data = data72, cutpoint = 0, bw = 30)
summary(rd4b)

low4b <- rd4b$ci[1,1]
up4b <- rd4b$ci[1,2]
coef4b <- rd4b$est[1]
se4b <- rd4b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ four | region, data = data72, cutpoint = 0, bw = 30))

# with covariates
cov4b <- RDestimate(factor ~ four | cntry + education + gndr + agea + foreignparent, data = data72, cutpoint = 0, bw = 30)
summary(cov4b)


# 21.09.2016
data81$five <- NA
data81$five <- as.numeric(data81$date) - 17065
hist(data81$five, breaks = 1000)
rd5 <- RDestimate(factor ~ five | cntry, data = data81, cutpoint = 0, bw = 30)
summary(rd5)

low5 <- rd5$ci[1,1]
up5 <- rd5$ci[1,2]
coef5 <- rd5$est[1]
se5 <- rd5$se[1]

# with region fixed effects
#summary(RDestimate(factor ~ five | region, data = data81, cutpoint = 0, bw = 30)) # NaNs produced

# with covariates
cov5 <- RDestimate(factor ~ five | cntry + education + gndr + agea + foreignparent, data = data81, cutpoint = 0, bw = 30)
summary(cov5)

data82$five <- NA
data82$five <- as.numeric(data82$date) - 17065
hist(data82$five, breaks = 1000)
rd5b <- RDestimate(factor ~ five | cntry, data = data82, cutpoint = 0, bw = 30)
summary(rd5b)

low5b <- rd5b$ci[1,1]
up5b <- rd5b$ci[1,2]
coef5b <- rd5b$est[1]
se5b <- rd5b$se[1]

# with region fixed effects
#summary(RDestimate(factor ~ five | region, data = data82, cutpoint = 0, bw = 30)) # NaNs produced

# with covariates
cov5b <- RDestimate(factor ~ five | cntry + education + gndr + agea + foreignparent, data = data82, cutpoint = 0, bw = 30)
summary(cov5b)

# 02.11.2016
data81$six <- as.numeric(data81$date) - 17107
hist(data81$six, breaks = 1000)
rd6 <- RDestimate(factor ~ six | cntry, data = data81, cutpoint = 0, bw = 30)
summary(rd6)

low6 <- rd6$ci[1,1]
up6 <- rd6$ci[1,2]
coef6 <- rd6$est[1]
se6 <- rd6$se[1]

# with region fixed effects
summary(RDestimate(factor ~ six | region, data = data81, cutpoint = 0, bw = 30))

# with covariates
cov6 <- RDestimate(factor ~ six | cntry + education + gndr + agea + foreignparent, data = data81, cutpoint = 0, bw = 30)
summary(cov6)

data82$six <- as.numeric(data82$date) - 17107
hist(data82$six, breaks = 1000)
rd6b <- RDestimate(factor ~ six | cntry, data = data82, cutpoint = 0, bw = 30)
summary(rd6b)

low6b <- rd6b$ci[1,1]
up6b <- rd6b$ci[1,2]
coef6b <- rd6b$est[1]
se6b <- rd6b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ six | region, data = data82, cutpoint = 0, bw = 30))

# with covariates
cov6b <- RDestimate(factor ~ six | cntry + education + gndr + agea + foreignparent, data = data82, cutpoint = 0, bw = 30)
summary(cov6b)

# 14.11.2016
data81$seven <- as.numeric(data81$date) - 17119
hist(data81$seven, breaks = 1000)
rd7 <- RDestimate(factor ~ seven | cntry, data = data81, cutpoint = 0, bw = 30)
summary(rd7)

low7 <- rd7$ci[1,1]
up7 <- rd7$ci[1,2]
coef7 <- rd7$est[1]
se7 <- rd7$se[1]

# with region fixed effects
summary(RDestimate(factor ~ seven | region, data = data81, cutpoint = 0, bw = 30))

# with covariates
cov7 <- RDestimate(factor ~ seven | cntry + education + gndr + agea + foreignparent, data = data81, cutpoint = 0, bw = 30)
summary(cov7)

data82$seven <- as.numeric(data82$date) - 17119
hist(data82$seven, breaks = 1000)
rd7b <- RDestimate(factor ~ seven | cntry, data = data82, cutpoint = 0, bw = 30)
summary(rd7b)

low7b <- rd7b$ci[1,1]
up7b <- rd7b$ci[1,2]
coef7b <- rd7b$est[1]
se7b <- rd7b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ seven | region, data = data82, cutpoint = 0, bw = 30))

# with covariates
cov7b <- RDestimate(factor ~ seven | cntry + education + gndr + agea + foreignparent, data = data82, cutpoint = 0, bw = 30)
summary(cov7b)

# 14.01.2017
data81$eight <- as.numeric(data81$date) - 17180
hist(data81$eight, breaks = 1000)
rd8 <- RDestimate(factor ~ eight | cntry, data = data81, cutpoint = 0, bw = 30)
summary(rd8)

low8 <- rd8$ci[1,1]
up8 <- rd8$ci[1,2]
coef8 <- rd8$est[1]
se8 <- rd8$se[1]

# with region fixed effects
summary(RDestimate(factor ~ eight | region, data = data81, cutpoint = 0, bw = 30))

# with covariates
cov8 <- RDestimate(factor ~ eight | cntry + education + gndr + agea + foreignparent, data = data81, cutpoint = 0, bw = 30)
summary(cov8)

data82$eight <- as.numeric(data82$date) - 17180
hist(data82$eight, breaks = 1000)
rd8b <- RDestimate(factor ~ eight | cntry, data = data82, cutpoint = 0, bw = 30)
summary(rd8b)

low8b <- rd8b$ci[1,1]
up8b <- rd8b$ci[1,2]
coef8b <- rd8b$est[1]
se8b <- rd8b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ eight | region, data = data82, cutpoint = 0, bw = 30))

# with covariates
cov8b <- RDestimate(factor ~ eight | cntry + education + gndr + agea + foreignparent, data = data82, cutpoint = 0, bw = 30)
summary(cov8b)


# 20.02.2017
data81$nine <- as.numeric(data81$date) - 17217
hist(data81$nine, breaks = 1000)
rd9 <- RDestimate(factor ~ nine | cntry, data = data81, cutpoint = 0, bw = 30)
summary(rd9)

low9 <- rd9$ci[1,1]
up9 <- rd9$ci[1,2]
coef9 <- rd9$est[1]
se9 <- rd9$se[1]

# with region fixed effects
summary(RDestimate(factor ~ nine | region, data = data81, cutpoint = 0, bw = 30))

# with covariates
cov9 <- RDestimate(factor ~ nine | cntry + education + gndr + agea + foreignparent, data = data81, cutpoint = 0, bw = 30)
summary(cov9)

data82$nine <- as.numeric(data82$date) - 17217
hist(data82$nine, breaks = 1000)
rd9b <- RDestimate(factor ~ nine | cntry, data = data82, cutpoint = 0, bw = 30)
summary(rd9b)

low9b <- rd9b$ci[1,1]
up9b <- rd9b$ci[1,2]
coef9b <- rd9b$est[1]
se9b <- rd9b$se[1]

# with region fixed effects
#summary(RDestimate(factor ~ nine | region, data = data82, cutpoint = 0, bw = 30)) # NaNs produced

# with covariates
cov9b <- RDestimate(factor ~ nine | cntry + education + gndr + agea + foreignparent, data = data82, cutpoint = 0, bw = 30)
summary(cov9b)

# 23.03.2017
data81$ten <- as.numeric(data81$date) - 17248
hist(data81$ten, breaks = 1000)
rd10 <- RDestimate(factor ~ ten | cntry, data = data81, cutpoint = 0, bw = 30)
summary(rd10)

low10 <- rd10$ci[1,1]
up10 <- rd10$ci[1,2]
coef10 <- rd10$est[1]
se10 <- rd10$se[1]

# with region fixed effects
summary(RDestimate(factor ~ ten | region, data = data81, cutpoint = 0, bw = 30))

# with covariates
cov10 <- RDestimate(factor ~ ten | cntry + education + gndr + agea + foreignparent, data = data81, cutpoint = 0, bw = 30)
summary(cov10)

data82$ten <- as.numeric(data82$date) - 17248
hist(data82$ten, breaks = 1000)
rd10b <- RDestimate(factor ~ ten | cntry, data = data82, cutpoint = 0, bw = 30)
summary(rd10b)

low10b <- rd10b$ci[1,1]
up10b <- rd10b$ci[1,2]
coef10b <- rd10b$est[1]
se10b <- rd10b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ ten | region, data = data82, cutpoint = 0, bw = 30))

# with covariates
cov10b <- RDestimate(factor ~ ten | cntry + education + gndr + agea + foreignparent, data = data82, cutpoint = 0, bw = 30)
summary(cov10b)

# 16.04.2017
data81$eleven <- as.numeric(data81$date) - 17272
hist(data81$eleven, breaks = 1000)
rd11 <- RDestimate(factor ~ eleven | cntry, data = data81, cutpoint = 0, bw = 30)
summary(rd11)

low11 <- rd11$ci[1,1]
up11 <- rd11$ci[1,2]
coef11 <- rd11$est[1]
se11 <- rd11$se[1]

# with region fixed effects
summary(RDestimate(factor ~ eleven | region, data = data81, cutpoint = 0, bw = 30))

# with covariates
cov11 <- RDestimate(factor ~ eleven | cntry + education + gndr + agea + foreignparent, data = data81, cutpoint = 0, bw = 30)
summary(cov11)

data82$eleven <- as.numeric(data82$date) - 17272
hist(data82$eleven, breaks = 1000)
rd11b <- RDestimate(factor ~ eleven | cntry, data = data82, cutpoint = 0, bw = 30)
summary(rd11b)

low11b <- rd11b$ci[1,1]
up11b <- rd11b$ci[1,2]
coef11b <- rd11b$est[1]
se11b <- rd11b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ eleven | region, data = data82, cutpoint = 0, bw = 30)) 

# with covariates
cov11b <- RDestimate(factor ~ eleven | cntry + education + gndr + agea + foreignparent, data = data82, cutpoint = 0, bw = 30)
summary(cov11b)


# 17.06.2017
data81$twelve <- as.numeric(data81$date) - 17334
hist(data81$twelve, breaks = 1000)
rd12 <- RDestimate(factor ~ twelve | cntry, data = data81, cutpoint = 0, bw = 30)
summary(rd12)

low12 <- rd12$ci[1,1]
up12 <- rd12$ci[1,2]
coef12 <- rd12$est[1]
se12 <- rd12$se[1]

# with region fixed effects
#summary(RDestimate(factor ~ twelve | region, data = data81, cutpoint = 0, bw = 30)) NaNs produced

# with covariates
cov12 <- RDestimate(factor ~ twelve | cntry + education + gndr + agea + foreignparent, data = data81, cutpoint = 0, bw = 30)
summary(cov12)

data82$twelve <- as.numeric(data82$date) - 17334
hist(data82$twelve, breaks = 1000)
rd12b <- RDestimate(factor ~ twelve | cntry, data = data82, cutpoint = 0, bw = 30)
summary(rd12b)

low12b <- rd12b$ci[1,1]
up12b <- rd12b$ci[1,2]
coef12b <- rd12b$est[1]
se12b <- rd12b$se[1]

# with region fixed effects
#summary(RDestimate(factor ~ twelve | region, data = data82, cutpoint = 0, bw = 30)) NaNs produced

# with covariates
cov12b <- RDestimate(factor ~ twelve | cntry + education + gndr + agea + foreignparent, data = data82, cutpoint = 0, bw = 30)
summary(cov12b)


# 18.01.2019
data91$thirteen <- as.numeric(data91$date) - 17914
hist(data91$thirteen, breaks = 1000)
rd13 <- RDestimate(factor ~ thirteen | cntry, data = data91, cutpoint = 0, bw = 30)
summary(rd13)

low13 <- rd13$ci[1,1]
up13 <- rd13$ci[1,2]
coef13 <- rd13$est[1]
se13 <- rd13$se[1]

# with region fixed effects
summary(RDestimate(factor ~ thirteen | region, data = data91, cutpoint = 0, bw = 30))

# with covariates
cov13 <- RDestimate(factor ~ thirteen | cntry + education + gndr + agea + foreignparent, data = data91, cutpoint = 0, bw = 30)
summary(cov13)

data92$thirteen <- as.numeric(data92$date) - 17914
hist(data92$thirteen, breaks = 1000)
rd13b <- RDestimate(factor ~ thirteen | cntry, data = data92, cutpoint = 0, bw = 30)
summary(rd13b)

low13b <- rd13b$ci[1,1]
up13b <- rd13b$ci[1,2]
coef13b <- rd13b$est[1]
se13b <- rd13b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ thirteen | region, data = data92, cutpoint = 0, bw = 30))

# with covariates
cov13b <- RDestimate(factor ~ thirteen | cntry + education + gndr + agea + foreignparent, data = data92, cutpoint = 0, bw = 30)
summary(cov13b)

### 18.01.2019 ITALY ONLY ###
data91$thirteen <- as.numeric(data91$date) - 17914
hist(data91$thirteen, breaks = 1000)

data91IT <- data91 %>% 
  filter(cntry == "IT")
rd13IT <- RDestimate(factor ~ thirteen, data = data91IT, cutpoint = 0, bw = 30)
summary(rd13IT)

low13IT <- rd13IT$ci[1,1]
up13IT <- rd13IT$ci[1,2]
coef13IT <- rd13IT$est[1]
se13IT <- rd13IT$se[1]

# with region fixed effects
summary(RDestimate(factor ~ thirteen | region, data = data91IT, cutpoint = 0, bw = 30))

# with covariates
cov13IT <- RDestimate(factor ~ thirteen | education + gndr + agea + foreignparent, data = data91IT, cutpoint = 0, bw = 30)
summary(cov13IT)

data92$thirteen <- as.numeric(data92$date) - 17914
hist(data92$thirteen, breaks = 1000)

data92IT <- data92 %>% 
  filter(cntry == "IT")
rd13bIT <- RDestimate(factor ~ thirteen, data = data92IT, cutpoint = 0, bw = 30)
summary(rd13bIT)

low13bIT <- rd13bIT$ci[1,1]
up13bIT <- rd13bIT$ci[1,2]
coef13bIT <- rd13bIT$est[1]
se13bIT <- rd13bIT$se[1]

# with region fixed effects
summary(RDestimate(factor ~ thirteen | region, data = data92IT, cutpoint = 0, bw = 30))

# with covariates
cov13bIT <- RDestimate(factor ~ thirteen | education + gndr + agea + foreignparent, data = data92IT, cutpoint = 0, bw = 30)
summary(cov13bIT)


# 04.12.2019
data91$fourteen <- as.numeric(data91$date) - 18234
hist(data91$fourteen, breaks = 1000)
rd14 <- RDestimate(factor ~ fourteen | cntry, data = data91, cutpoint = 0, bw = 30)
summary(rd14)

low14 <- rd14$ci[1,1]
up14 <- rd14$ci[1,2]
coef14 <- rd14$est[1]
se14 <- rd14$se[1]

# with region fixed effects
summary(RDestimate(factor ~ fourteen | region, data = data91, cutpoint = 0, bw = 30))

# with covariates
cov14 <- RDestimate(factor ~ fourteen | cntry + education + gndr + agea + foreignparent, data = data91, cutpoint = 0, bw = 30)
summary(cov14)

data92$fourteen <- as.numeric(data92$date) - 18234
hist(data92$fourteen, breaks = 1000)
rd14b <- RDestimate(factor ~ fourteen | cntry, data = data92, cutpoint = 0, bw = 30)
summary(rd14b)

low14b <- rd14b$ci[1,1]
up14b <- rd14b$ci[1,2]
coef14b <- rd14b$est[1]
se14b <- rd14b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ fourteen | region, data = data92, cutpoint = 0, bw = 30))

# with covariates
cov14b <- RDestimate(factor ~ fourteen | cntry + education + gndr + agea + foreignparent, data = data92, cutpoint = 0, bw = 30)
summary(cov14b)


# 17.12.2021

data101$fifteen <- as.numeric(data101$date) - 18978
hist(data101$fifteen, breaks = 1000)
rd15 <- RDestimate(factor ~ fifteen | cntry, data = data101, cutpoint = 0, bw = 30)
summary(rd15)

low15 <- rd15$ci[1,1]
up15 <- rd15$ci[1,2]
coef15 <- rd15$est[1]
se15 <- rd15$se[1]

# with region fixed effects
summary(RDestimate(factor ~ fifteen | region, data = data101, cutpoint = 0, bw = 30))

# with covariates
cov15 <- RDestimate(factor ~ fifteen | cntry + education + gndr + agea + foreignparent, data = data101, cutpoint = 0, bw = 30)
summary(cov15)

data102$fifteen <- as.numeric(data102$date) - 18978
hist(data102$fifteen, breaks = 1000)
rd15b <- RDestimate(factor ~ fifteen | cntry, data = data102, cutpoint = 0, bw = 30)
summary(rd15b)

low15b <- rd15b$ci[1,1]
up15b <- rd15b$ci[1,2]
coef15b <- rd15b$est[1]
se15b <- rd15b$se[1]

# with region fixed effects
summary(RDestimate(factor ~ fifteen | region, data = data102, cutpoint = 0, bw = 30))

# with covariates
cov15b <- RDestimate(factor ~ fifteen | cntry + education + gndr + agea + foreignparent, data = data102, cutpoint = 0, bw = 30)
summary(cov15b)


### 17.12.2021 ITALY ONLY ###

data101$fifteen <- as.numeric(data101$date) - 18978
hist(data101$fifteen, breaks = 1000)

data101IT <- data101 %>% 
  filter(cntry == "IT")
rd15IT <- RDestimate(factor ~ fifteen, data = data101IT, cutpoint = 0, bw = 30)
summary(rd15IT)

low15IT <- rd15IT$ci[1,1]
up15IT <- rd15IT$ci[1,2]
coef15IT <- rd15IT$est[1]
se15IT <- rd15IT$se[1]

# with region fixed effects
summary(RDestimate(factor ~ fifteen | region, data = data101IT, cutpoint = 0, bw = 30))

# with covariates
cov15IT <- RDestimate(factor ~ fifteen | education + gndr + agea + foreignparent, data = data101IT, cutpoint = 0, bw = 30)
summary(cov15IT)


data102$fifteen <- as.numeric(data102$date) - 18978
hist(data102$fifteen, breaks = 1000)

data102IT <- data102 %>% 
  filter(cntry == "IT")
rd15bIT <- RDestimate(factor ~ fifteen, data = data102IT, cutpoint = 0, bw = 30)
summary(rd15bIT)

low15bIT <- rd15bIT$ci[1,1]
up15bIT <- rd15bIT$ci[1,2]
coef15bIT <- rd15bIT$est[1]
se15bIT <- rd15bIT$se[1]

# with region fixed effects
summary(RDestimate(factor ~ fifteen | region, data = data102IT, cutpoint = 0, bw = 30))

# with covariates
cov15bIT <- RDestimate(factor ~ fifteen | education + gndr + agea + foreignparent, data = data102IT, cutpoint = 0, bw = 30)
summary(cov15bIT)

sink() 
